Breast cancer
Mariana Pasqualini
10/15/2020
Motivação do classificador de câncer de mama
O objetivo desse notebook é aplicar técnicas de classificação aos dados do dataset Breast cancer Wisconsin, disponível aqui. O interesse nesse dataset vem do outubro rosa, mês de conscientização para informar, alertar e previnir o câncer de mama. De acordo com o Instituto Nacional de Câncer (INCA), em 2019 foram estimados 59700 novos casos da doença no Brasil.
Quanto antes, melhor!
Esse estudo é apenas para fins didáticos e de aprendizado, que não necessariamente reflete a realidade ou uma opinião médica.
Pacotes
library(tidyverse)
library(janitor)
library(GGally)
library(caret)
library(printr)
library(corrplot)
library(factoextra)
library(pROC)
set.seed(889)
source("analysis.r")
## [1] "AUC on Test 0.999689440993789"Explorando os dados: Correlações e separabilidade das classes
Mean
É possível ver rapidamente que algumas variáveis estão correlacionadas fortemente entre si.
Algumas variáveis apresentam correlação positiva quase perfeita (muito próximo de 1), como
radius_mean e perimeter_mean (0.998)
radius_mean e area_mean (0.987)
perimeter_area e area_mean (0.987)
Outras variáveis também apresentam uma correlação forte (acima de 0.5), como
concavity_mean e compactness_mean (0.883)
concave.points_mean e concavity_mean (0.921)
concave.points_mean e perimeter_mean (0.851)
Quanto à separabilidade dos dados, observando os histogramas, parece que as variáveis radius_mean, perimeter_mean, area_mean e concavity_mean ajudam a melhor identificar as populações benigno e maligno, mas ainda assim vemos sobreposições e nenhuma variável parece separar perfeitamente essas classes.
Standard Error
As correlações do erro padrão dessas variáveis são, de maneira geral, um pouco mais baixas, mas ainda sim temos correlações fortes nas mesmas variáveis da média.
Quanto à separabilidade das classes, as variáveis parecem não ajudar a identificar as duas populações.
Worst
Os scatterplots dessas variáveis, que são as médias dos maiores valores, também são bem parecidas com as observadas para as médias.
Por meio das variáveis concave.points_worst, radius_worst e perimeter_worst parecem ser as que mais ajudam separar as classes do diagnóstico.
Quando uma variável tem correlação perfeita, sabendo o valor da variável X, conseguimos prever o de Y. Com isso, podemos escrever X em função de Y (ou vice-versa). E observamos nas matrizes acima que algumas variáveis têm correlação bem alta.
Trabalhar com dados altamente correlacionados:
- Não acrescenta informação da variabilidade dos dados
- Os modelos podem levar mais tempo que o necessário para serem ajustados
Para ajudar nesse problema, vamos recorrer à técnicas de redução de dimensionalidade.
PCA: Reduzindo a dimensão
Análise de componentes principais é uma técnica de fatorização de matriz em que é possível explicar (parte da) a variabilidade dos dados por meio de combinações lineares não correlacionadas dos dados originais. De maneira geral, a ideia de fatorar uma matriz é escrever uma matriz “complicada” no produto de duas matrizes “mais simples”.
Vamos definir as componentes principais por meio de operações matriciais, para conhecer todo o processo.
- Estimando a matriz de correlação:
O primeiro passo é estimar a matriz de correlação dos nossos dados. Na verdade, nessa etapa usa-se a matriz de covariâncias, mas um problema que pode surgir é que algumas variância são maiores que outras, e em casos muito discrepantes isso distorce as componentes principais. Para evitar isso, o procedimento é: padronizar a matriz de covariâncias e estimar as componentes a partir dessa matriz padronizada. Mas isso é o mesmo que usar a matriz de correlação e por isso vou seguir por ela.
upper <- cor_matrix
upper[upper.tri(cor_matrix)] <- NA
knitr::kable(upper[1:5, 1:5], caption = "Matriz de correlação: 5 primeiras linhas/colunas")| radius_mean | texture_mean | perimeter_mean | area_mean | smoothness_mean | |
|---|---|---|---|---|---|
| radius_mean | 1.0000000 | NA | NA | NA | NA |
| texture_mean | 0.3237819 | 1.0000000 | NA | NA | NA |
| perimeter_mean | 0.9978553 | 0.3295331 | 1.0000000 | NA | NA |
| area_mean | 0.9873572 | 0.3210857 | 0.9865068 | 1.0000000 | NA |
| smoothness_mean | 0.1705812 | -0.0233885 | 0.2072782 | 0.1770284 | 1 |
- Obtendo os autovalores:
eig$values
## [1] 1.328161e+01 5.691355e+00 2.817949e+00 1.980640e+00 1.648731e+00
## [6] 1.207357e+00 6.752201e-01 4.766171e-01 4.168948e-01 3.506935e-01
## [11] 2.939157e-01 2.611614e-01 2.413575e-01 1.570097e-01 9.413497e-02
## [16] 7.986280e-02 5.939904e-02 5.261878e-02 4.947759e-02 3.115940e-02
## [21] 2.997289e-02 2.743940e-02 2.434084e-02 1.805501e-02 1.548127e-02
## [26] 8.177640e-03 6.900464e-03 1.589338e-03 7.488031e-04 1.330448e-04
plot(eig$values)
abline(h = 1)
Para estimar o número de componentes principais ideal, um critério prático é observar os autovalores maiores ou iguais a 1, pois dessa forma as combinações lineares explicam a variância de uma variável original (padronizada) do dataset. Observando o scree-plot acima, é k = 6 é o número de componentes principais adequados para explicar a estrutura da variância dos dados.
- Obtendo os autovetores e escrevendo a 1ª componente principal:
Observando os autovetores associados aos seis autovalores identificados (apenas as 6 primeiras linhas de 30)
| V1 | V2 | V3 | V4 | V5 | V6 |
|---|---|---|---|---|---|
| -0.2189024 | -0.2338571 | -0.0085312 | 0.0414090 | 0.0377864 | 0.0187408 |
| -0.1037246 | -0.0597061 | 0.0645499 | -0.6030500 | -0.0494689 | -0.0321788 |
| -0.2275373 | -0.2151814 | -0.0093142 | 0.0419831 | 0.0373747 | 0.0173084 |
| -0.2209950 | -0.2310767 | 0.0286995 | 0.0534338 | 0.0103313 | -0.0018877 |
| -0.1425897 | 0.1861130 | -0.1042919 | 0.1593828 | -0.3650885 | -0.2863745 |
| -0.2392854 | 0.1518916 | -0.0740916 | 0.0317946 | 0.0117040 | -0.0141309 |
Olhando um pouco para a teoria, a j-ésima componente principal é dada por:
\[ Ŷ_j = ê_{j1}X_{1} + ê_{j2}X_{2} + ... + ê_{jp}X_{p} \]
E podemos escrever, para os dados do dataset em questão, a primeira componente principal da seguinte forma:
\[ Ŷ_1 = -0.2189\text{radius_mean} -0.1037\text{texture_mean} - 0.227\text{perimeter_mean} + ... - 0.132\text{fractal_dimension_worst} \]
E assim em diante…
Mas como visualizações são mais interessantes, vamos construir gráficos!
Visualizando as componentes principais
Bom, como a função fviz_pca_biplot do pacote factoextra faz um gráfico apropriado a partir de um objeto da classe PCA, vou utilizar a função princomp (do stats, que faz parte do R base) para criar esse objeto. Essa função foi escolhida porque ela usa o teorema da decomposição espectral, que é a abordagem teórica utilizada na descrição acima. Caso a abordagem fosse decomposição em valores singulares (SVD), a função apropriada seria a prcomp (também do stats).
pca <- princomp(data_processed[,-1], cor = TRUE, scores = TRUE)
fviz_pca_biplot(pca, col.ind = data_processed$diagnosis, col="black",
palette = "Dark2", geom = "point", repel=TRUE,
legend.title="Diagnóstico", addEllipses = TRUE) +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) 
A primeira componente explica 44.3% da variância dos dados. É possível observar que existem variáveis correlacionadas negativamente, como smoothness_se e radius_mean, justamente os vetores mais distantes entre si. Também vemos que há ali um agrupamento das classes Benigno/Maligno.
Como o PCA foi utilizado para reduzir a dimensionalidade dos dados, os scores das 6 primeiras dimensões serão utilizados para ajustar o modelo de classificação.
Classificando
Organizando os dados que serão utilizados no modelo:
data_scores <-
bind_cols(data_processed$diagnosis, as.data.frame(pca$scores[,1:6])) %>%
clean_names() %>%
rename(diagnosis = x1)## New names:
## * NA -> ...1Dividindo os dados em conjuntos de treino e teste (70/30):
train <- data_scores %>% sample_frac(.7)
test <- data_scores %>% anti_join(train)Modelo 01: Regressão logística
Ajustando o modelo
logistic_model1 <- glm(diagnosis ~ ., family = binomial, data = train)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredsummary(logistic_model1)##
## Call:
## glm(formula = diagnosis ~ ., family = binomial, data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6845 -0.0496 -0.0045 0.0012 3.4221
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.2843 0.3626 -0.784 0.43300
## comp_1 2.9268 0.5107 5.731 9.96e-09 ***
## comp_2 1.6343 0.3493 4.679 2.89e-06 ***
## comp_3 0.4722 0.2297 2.056 0.03982 *
## comp_4 -0.7846 0.2670 -2.939 0.00329 **
## comp_5 -1.3397 0.4449 -3.011 0.00260 **
## comp_6 -0.3746 0.2821 -1.328 0.18411
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 533.015 on 397 degrees of freedom
## Residual deviance: 57.436 on 391 degrees of freedom
## AIC: 71.436
##
## Number of Fisher Scoring iterations: 9A variável comp_6 parece não contribuir significativamente para classificarmos o diagnóstico (utilizando \(\alpha = 0.05\)). Um novo modelo pode ser ajustado removendo essa variável, já que variáveis irrelevantes podem deteriorar a taxa de erro do modelo.
logistic_model2 <- glm(diagnosis ~ comp_1 + comp_2 + comp_3 + comp_4 + comp_5,
family = binomial, data = train)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredsummary(logistic_model2)##
## Call:
## glm(formula = diagnosis ~ comp_1 + comp_2 + comp_3 + comp_4 +
## comp_5, family = binomial, data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7274 -0.0542 -0.0060 0.0016 3.4014
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.2408 0.3643 -0.661 0.50861
## comp_1 2.8145 0.4734 5.945 2.77e-09 ***
## comp_2 1.5617 0.3314 4.713 2.44e-06 ***
## comp_3 0.3441 0.1971 1.746 0.08086 .
## comp_4 -0.7896 0.2631 -3.002 0.00269 **
## comp_5 -1.2132 0.3975 -3.052 0.00227 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 533.015 on 397 degrees of freedom
## Residual deviance: 59.161 on 392 degrees of freedom
## AIC: 71.161
##
## Number of Fisher Scoring iterations: 9E como fica no conjunto de teste?
Queremos prever novas observações, então vamos lá:
pred <-
predict(logistic_model2, test, type = "response") %>%
as.data.frame() %>%
rename(., "prob" = ".") %>%
bind_cols(test) %>%
mutate(predicted = as.factor(ifelse(prob > 0.5, "Malignant", "Benign")))As probabilidades são dadas na forma \(P(Y = 1 | X)\), e olhando para a dummy definida como 1
contrasts(data_scores$diagnosis)| Malignant | |
|---|---|
| Benign | 0 |
| Malignant | 1 |
olhamos para os valores como sendo a probabilidade do tumor ser maligno.
Avaliando a acurácia do modelo
cmat <- conf_mat(table(pred$predicted, pred$diagnosis))
autoplot(cmat, type = "heatmap")
pred %>% metrics(diagnosis, predicted)| .metric | .estimator | .estimate |
|---|---|---|
| accuracy | binary | 0.9766082 |
| kap | binary | 0.9473765 |
A acurácia do modelo é de 99.4% e o coeficiente de Kappa 98.7%, que nos indica uma concordância entre o predito e os valores verdadeiros. É importante lembrar que a acurácia é uma métrica que deve ser usada com cuidado, pois ela pode não representar bem a situação.
Na área da saúde, as medidas mais utilizadas são a especificidade e a sensibilidade, comuns em testes diagnósticos. Em termos práticos, a sensibilidade avalia a capacidade do classificador detectar o tumor maligno quando ele está presente. Já a especificidade nos dá a probabilidade de um teste dar negativo (no caso, benigno) quando, de fato, não há doença. A sensibilidade do classificador é de 98,2% e a especifidade 100%.
Quando estamos em um cenário de classificação, uma forma de checarmos a taxa de erro associada ao conjunto de teste é estimada da seguinte forma:
\(Media(I(y_{0} \neq \hat{y_0}))\)
Ou seja, é a média de uma variável indicadora que assume 1 quando a classe predita é igual a classe verdadeira e 0, caso contrário. Queremos o menor valor possível para essa média. Assim, considerando o problema, temos:
mean(pred$diagnosis != pred$predicted)## [1] 0.02339181Ou seja, o modelo de regressão logística previu incorretamente, em média, apenas 0.6% das vezes.
ROC e AUC:
predr <- prediction(pred$prob, pred$diagnosis)
perf <- performance(predr, "tpr", "fpr")
plot(perf)
auc <- performance(predr, measure = 'auc')
auc <- unlist(slot(auc, 'y.values'))
print(paste('AUC on Test', auc))
## [1] "AUC on Test 0.996428571428571"